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ABSTRACT 


A  technique  for  reducing  the  number  of  iterations  necessary 
for  solving  linear  programs  using  the  primal-dual  algorithm  is 
presented.   It  appears  that  the  new  method  will  also  decrease  the 
number  of  iterations  over  any  other  simplex  algorithm.   A  FORTRAN 
program  incorporating  the  technique,  as  well  as  some  comparative 
computational  results  are  given. 
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SECTION  I 
INTRODUCTION 

The  primal-dual  algorithm,  as  developed  by  Ford  and  Fulkerson, 
theoretically  permitted  faster  convergence  in  solving  linear  programs 
over  such  methods  as  the  two-phase  or  Charnes '  M-method.   Instead  of 
first  driving  the  artificial  variables  to  zero  to  obtain  a  feasible 
solution  and  then  working  to  optimality  in  a  second  phase,  the  primal- 
dual  works  towards  optimality  and  feasibility  in  the  primal  and  dual 
problems  at  the  same  time.   This  is  accomplished  during  what  would 
be  the  first  phase  for  the  other  procedures.   However,  in  actual 
computational  experience,  the  primal-dual  algorithm  has  not  achieved 
faster  convergence  over  the  other  simplex  procedures. 

This  paper  presents  a  modification  of  the  primal-dual  algorithm 
that  does  accomplish  faster  convergence  than  the  usual  primal-dual 
algorithm  and  perhaps  also  the  other  simplex  methods. 


Professor  Harold  Greenberg  provided  the  idea  for  the  new 
technique. 


SECTION  II 


THE  STANDARD  PRIMAL- DUAL  ALGORITHM 


Using  matrix  notation,  the  linear  programming  problem  is: 

Minimize  z  =  z  +  CX 
o 

Given    AX  =  H 

x.^0     j=l, . . . ,n 
J 

where  z   is  a  constant  and 
o 

C  is  a  1  x  n  vector  with  components  c.  ^  0,  j  =  l,...,n 

H  is  an  m  x  1  vector  with  components  h.  ^  0,  i  =  l,...,m 

A  is  an  m  x  n  matrix  with  elements  a.  .   i  =  l,...,m 

l  >  J  > 

j  =  1,  —  ,n 

X  is  an  n  x  1  vector  with  components  x.,  j  =  l,...,n 

2 

The  standard  primal-dual,  as  it  appears  in  Dantzig,   consists 

in  augmenting  the  constraint  equations  with  artificial  variables  as, 


AX  +  Y  =  H  (1) 

yi  ^  0         i  =  1, . .  .,m 


where  Y  is  an  m  x  1  vector  with  components  consisting  of  the  arti- 
ficial variables  y.;  and  forming 


w  =  w  +  DX 
o 

where 

m 

w  =  S  h. 
o   .     1 
i=l 

and  D  is  a  1  x  n  vector  with  components  d.,  j  =  l,...,n  where 


2 
George  B.  Dantzig,  Linear  Programming  and  Extensions 

(Princeton:   Princeton  University  Press,  1963),  pp.  247-53. 


m 
d.  =  -  S   a. 
J     i=l  1>1 

It  is  noted  that  both  the  c.  and  the  h.  must  be  positive 

J  i 

prior  to  starting  the  computations.   Negative  h.  may  be  made 
positive  by  multiplying  the  associated  constraint  equations  by 
a  minus  one.   The  method  for  correcting  negative  c.  is  discussed 
in  the  "Description  of  the  Computer  Program"  section  of  this 
paper. 

The  problem  is  to  find  x.  ^  0,  w  =  0,  that  minimizes  z 
subject  to  (1).  At  any  stage  of  the  iterative  procedure  the 
equations  have  the  form: 

XB  +  AtX  =  Ht 

Y  +  A£X  =  H2  (2) 

DX  =  w  -  w 
o 

CX  =  z  -  z 
o 

where  X   is  a  vector  representing  the  current  basic  variables 

D 

and  X  is  a  vector  representing  the  current  non-basic  variables. 

The  numerical  subscripts  1  and  2  represent  partitioned  matrices 

and  vectors.   Matrices  and  vectors  with  bars  are  the  corresponding 

original  matrices  and  vectors  after  the  usual  pivot  operations 

within  the  primal  and  dual  problems.  When  an  artificial  variable 

becomes  non-basic,  it  is  dropped  from  further  consideration. 

This  results  in  the  A  matrix  and  the  H  vector.   Each  iteration 

in  the  standard  primal-dual  algorithm  consists  of  the  following: 

a)   Finding  the  most  negative  d.  value  for  those  indices 

i  having  "c.  =0. 
J 


b)  If  there  are  no  indices  with  this  property,  the 

multiplier  k  =  min  [c./  -  d.]  is  used  to  obtain  at 

d  .<0   J     J 
J 
least  one  c.  equal  to  zero  by 

c.  =  c.  +  k(d.) 
J    J      J 

c)  The   index  r   =    i    is   obtained   for   the  minimum  d.   with 

J 

c.  =  0.   Then  x   is  selected  to  be  basic. 
J  r 

d)  Then  min  (h./a.   )  is  found  for  a.    >  0.   This  will 

l   l,r  l,r 

select  a  pivot  row  in  either  the  L  or  Y  sets  of 

o 

equations . 

e)  The  simplex  pivot  method  is  then  used.   If  the  pivot 
row  is  in  the  Y  set  of  equations  the  artificial  variable 
in  that  row  is  dropped.   Thus,  in  any  case,  one  of  the 
forms  in  (2)  is  achieved. 

When  w  =  0,  the  basic  solution  is  optimal  and  the  problem 
is  solved. 

The  actual  computations  should  include  rules  to  prevent 
cycling,  such  as  the  perturbation  technique  or  the  use  of 
lexicographic  ordering. 
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SECTION  III 

DESCRIPTION  OF  THE  NEW  TECHNIQUE 

The  primal-dual  algorithm  is  modified  in  a  two-step  procedure. 
In  the  first  step,  a  pivot  row  in  the  Y  set  of  equations  is  selected, 
By  only  considering  elimination  of  artificial  variables,  after  m 
or  less  iterations,  w  equals  zero.  At  this  point  the  equations 
have  the  form : 

Xro  +  AX  =  H 
B 

z   +  CX  =  z 
o 

where  all  c.  are  positive.   If  all  the  h.  are  positive,  the  solution 

is  optimal.  We  then  stop.   However,  in  most  cases,  this  stage  of 

the  problem  will  contain  one  or  more  negative  h..   In  this  case 

the  subscript  "it"  is  assigned  to  the  row  with  the  most  negative 

h..   The  following  procedure  is  then  used  as  the  second  step  of  the 

modification: 

a)  The  equation  containing  h.   is  multiplied  through  by  minus 

one.   This  produces  h.   >  0. 

it 

b)  The  "it"  equation  is  then  added  to  all  equations  with 
negative  h. .  All  basic  variables  except  the  "it"  equation 
now  become  feasible. 

c)  To  maintain  the  canonical  form  an  artificial  variable  is 

added  to  the  "it"  equation.  A  new  infeasibility  equation 

with  w  =  h. .  is  formed  from  the  "it"  equation, 
o    it  n 

The  process  now  reverts  to  the  primal-dual  algorithm  to  achieve 
w  equal  to  zero.   Since  the  canonical  form  has  been  maintained  and 
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all  c.  have  been  kept  positive  throughout  the  process,  once 
w  equals  zero,  the  solution  is  optimal. 
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SECTION  IV 

COMPUTATIONAL  RESULTS 

Several  problems  were  run  on  an  IBM  360/67  using  both  the 
standard  primal-dual  and  the  primal-dual  with  the  new  technique, 
The  problems  used  were  test  problems  from  the  SHARE  Linear 
Programming  Project  and  two  highly  degenerate  assignment  type 
problems.   Results  are  given  in  Table  I  with  the  two  assignment 
problems  listed  first. 

TABLE  I 
TABLE  OF  COMPARATIVE  RESULTS 


PRIMAL-DUAL  WITH 


SIZE 


STANDARD 
PRIMAL-DUAL 


SIMPLEX    THE  NEW  TECHNIQUE 
PROBLEM     (m  x  n)    ITERATIONS   Iterations   Seconds   Iterations   Seconds 


Assignment  74  x  150 

Assignment  34  x  376 

SHARE  1A  33  x  64  70 

SHARE  ID  27  x  45  61 

SHARE  IE  31  x  106  66 

SHARE  IF  66  x  135  184 

SHARE  2A  30  x  103  114 


89 
74 
39 
38 
96 
125 
79 


32 

12 

7 

30 
86 
29 


83 

97 

39 

38 

136 

173 

119 


29 
28 

7 

42 

114 

43 


In  the  two  assignment  problems  the  primal-dual  with  the  new  technique 
resulted  in  optima  different  from  that  of  the  standard  primal-dual  since 


Philip  Wolfe  and  Leola  Cutler,  "Experiments  in  Linear  Programming," 
Recent  Advances  in  Mathematical  Programming  (New  York:  McGraw-Hill 
Book  Co.,  Inc.,  1963),  p.  198. 
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these  problems  have  non-unique  solutions.   The  two  methods  resulted 
in  the  same  optima  in  all  other  problems. 

No  accurate  timing  data  is  available.   The  timing  data  in  the 
above  table  though  does  give  an  indication  of  the  reduction  in 
execution  times  possible  with  the  new  technique.   However,  the 
modified  method  does  perform  less  work  than  the  usual  method 
because  in  the  first  step  ratios  and  comparisons  are  only  taken 
for  the  artificial  rows.   The  restart  of  the  problem  in  the  second 
step  requires  only  an  insignificant  amount  of  time.   Thus,  in  the 
table,  the  number  of  iterations  represents  a  measure  of  improvement 
of  the  new  method.   The  new  technique  shows  significant  reductions 
in  the  number  of  iterations  from  the  two-phase  standard  simplex 
method  in  four  out  of  the  five  problems  where  a  comparison  is 
possible.  When  comparing  the  new  technique  with  the  standard 
primal-dual,  a  similar  reduction  is  evident  in  the  larger  problems. 
It  therefore  appears  that  the  new  method  presents  a  significant 
improvement  when  the  number  of  variables  is  large. 
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SECTION  V 
DESCRIPTION  OF  THE  COMPUTER  PROGRAM 

The  accompanying  program  (see  Appendix)  is  written  in  the 
FORTRAN  IV  language,  and  was  used  to  run  all  test  problems  with 
the  exception  of  the  two  assignment  type  problems  (since  these 
only  deal  with  zeros  or  ones  in  the  A  matrix,  considerable 
reduction  in  the  core  memory  requirement  could  be  made  with  a 
slight  modification  of  the  program).   The  program  will  handle  up 
to  69  constraint  equations  and  up  to  149  variables,  and  requires 
about  160,000  bits  of  core  storage.   Columns  may  be  designated 
by  alphanumeric  names,  but  rows  must  be  consecutive  integer 
values  not  to  exceed  the  value  69.   Only  the  non-zero  elements 
of  A,  C,  and  H  need  be  read  in.  All  h.  must  be  positive,  but 
negative  c  .  are  allowed. 

Names  used  in  this  program  conform  to  those  already  outlined 
in  this  paper.   Those  not  previously  discussed  are  defined  as: 

4 

FF  is  a  multiplicative  constant  used  to  insure  that  all  a.  . 

are  integer. 
EP  is  a  constant  used  in  the  round-off  error  control. 
DI  is  the  cumulative  multiplication  of  the  pivot  elements. 
CE  is  a  vector  identical  to  the  initial  C,  and  is  used  to 

calculate  z. 
AX  is  the  pivot  column  in  the  Revised  Simplex  tableau. 
CNAME  is  a  vector  of  column  names. 
DUAL  is  the  vector  of  dual  variables. 
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DB  is  a  vector  of  simplex  multipliers  used  to  calculate  the  d.. 
R  is  a  vector  used  to  record  whether  each  variable  is  basic 

or  artificial. 
E  is  a  computed  scalar  used  in  the  updating  and  round-off 

control. 
P  is  a  vector  of  multipliers  used  in  the  actual  pivot  operation. 
ITR  indicates  the  number  of  iterations  performed. 
INDEX  controls  the  use  of  the  new  technique  (INDEX  equals  zero 

yields  the  standard  primal-dual,  and  INDEX  equals  one 

incorporates  the  new  technique) . 
IB  is  a  vector  of  variable  designators  (a  negative  number 

indicates  an  artificial  variable,  and  a  positive  number 

indicates  a  basic  variable). 
All  other  names  are  defined  in  the  program. 

Once  the  data  has  been  read  in,  the  c.  are  tested  for  negative 
values.   If  all  c.  are  positive,  the  primal-dual  commences.   However, 
if  there  is  one  or  more  negative  c.  the  following  procedure  is  fol- 


lowed: 


a)  RMIN  is  set  equal  to  the  minimum  c. 

J 

b)  n  is  set  equal  to  n  +  1  and  m  is  set  equal  to  m  +  1. 

c)  The  coefficients  of  all  variables  in  the  added  constraint 
equation  indicated  in  b)  above  are  set  equal  to  FF  including 
x  .   (Variable  x  does  not  appear  in  any  other  constraint 
equation.) 

d)  Then  c   is  set  equal  to  minus  RMIN. 

n         n 
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e)   The  right  hand  side  value,  h  ,  is  entered  as: 

h  =  k(FF) 
m 

where  k  is  a  number  larger  than  the  sum  of  the  x. 

J 

in  the  optimal  tableau. 
The  problem  is  now  modified  so  it  may  be  handled  by  use  of  the 
primal-dual  algorithm. 

The  remainder  of  the  program  is  a  FORTRAN  representation  of 
the  steps  outlined  in  the  sections  "The  Standard  Primal-Dual 
Algorithm"  and  "Description  of  the  New  Technique."  However,  for 
computer  calculations,  an  essential  addition  to  the  procedures  is 
some  control  over  round-off  error. 

Round-off  error  develops  when  two  large  numbers,  say  y  and  z, 
(differing  only  slightly  from  each  other),  are  subtracted  from 
each  other  such  as: 

x  =  y  -  z 
In  this  program  the  absolute  value  of  x  divided  by  y  is  compared 
to  EP.   If  it  is  less  than  EP,  x  is  set  equal  to  zero.   As  a 
further  control,  the  values  of  DB,  B,  and  H  are  updated  each 
iteration. 
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C      PROGRAM  TO  SOLVE  PRIMAL-DUAL  LINEAR 

C      PROGRAMMING  PROBLEMS 

C 

•    IMPLICIT    REAL*8    (A-H,0-Z) 

COMMON  A (70, 150), 8(70.70) ,C( 150) ,D( 150) ,H(70) ,AX( 7C) 

COMMON  0B(7C) ,  DUAL ( 70 ) , R ( 1 50  I  , C E  ( 150),E,W,Z,DI,FF,EP 

COMMON  CNAME( 150) 

COMMON  18(70), M,N,JS,IR, INDEX, ITR,M1,N1 

KEAL*6  NINE/'9999      •/ 

FF=10000. 

EP=0.000C01 

DI  =  1. 

M=0 

DO  1  1=1,70 

H( I)=0. 

DO  1  J=1,15C 

1  A( I, J)=0. 

DO  2  1=1,150 
CNAMEd  )  =  C. 
C6(I  )=C. 

2  C(I)=0. 
C 

C      READING  IN  COEFICIENTS  OF  CONSTRAINT  EOUATIHNS 
C 

DO  5  N=l,1000 

READ  (5,3)  CB,  J,  R8 

IF(CB.EO.NINE)  GO  TO  16 

3  FORMAT  (6X,  A6,  I  6 ,  F10.4) 
DO  4  1=1,150 

IF(CNAME( I) .EO.CB)  GO  TO  30 
IF(CNAME(  I) .EO.C. )  GO  TO  30 
GO  TO  4 
30  A(J,I)=RB*FF 
CNAME( I )=CB 

WRITE  (6,102)  I,  J,  A(J,I) 
102  FORMAT  (IX,  2110,  FIR. 9) 
IF(J.LE.M)  GO  TO  5 
M=J 
GO  TO  5 

4  CONTINUE 

5  CONTINUE 
C 

C      READING  IN  RIGHT  HAND  SIDE  VALUES 
C 

16  DO  6  N=1,1C0 

READ  (5,14)  I,  RB 
14  FORMAT  ( 118,  F10.4) 
IFU.GT.999)  GO  TO  17 
H(I)=RB*FF 
WRITE  (6,14)  I,  H(I) 

6  CONTINUE 
C 

C      READING  IN  COSTS 


C 


17  DO  9  J=1,2CC 

READ  (5,7)  CC,  II,  RB 
7  FORMAT  (6X,  A6,  16,  F10.4) 

IF(CC.EQ.NINE)  GO  TO  18 

DO  11  K=l,150 

IF(CC.EQ.CNAME(K) )  GO  TO  12 

GO  TO  11 
12  C(K)=RB*FF 
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13 

11 

9 
18 


20 

21 

101 


103 
104 


CE(K)=RB 
WRITE  (6,13)  K, 
FORMAT  (6X,  216, 
GO  TO  9 
CONTINUE 
CONTINUE 
CONTINUE 
DO  20  1=1,150 
IF(CNAME( I) .GT 
N=I-1 
GO  TO  21 
CONTINUE 

WRITE  (6,101)  M,  N 
FORMAT  (10X,  'NUMBER 
l'NUMBER  OF  COLUMNS=», 
WRITE  (6,103)  (CNAME( 
FORMAT  (IX,  20A6) 
WRITE  (6,104)  (C(J), 
FORMAT  (IX,  10E12.4) 
INDEX=1 
N1=N+1 
M1=M+1 
CALL  LIP 
STOP 
END 


II,  C(K) 
F15.5) 


C.  )  GO  TO  20 


OF  ROWS=' 

15) 
I),  1=1,175) 

J=1,N) 


15,  15X 


21 


c 

c 


SUBROUTINE  LIP 


•  IMPLICIT  RFAL*8  <A-H,0-Z) 

COMMON  A(7C,15C),R(70,70) ,C ( 1 50 ) , D< 1 50 ) ,H(7C) ,AX( 70) 

COMMON  OR (70) ,DUAL(70),R(150) ,CF(150),F,W,Z,DI,FF,EP 

COMMON  CNAME( 150) 

COMMON  IB(70),M,N,JS, I R, INDEX, ITR,M1,N1 

E=0.5/DI 
C 

C      MODIFYING  THE  PROBLEM  IF  THE  OBJECTIVE 
C      FUNCTION  CONTAINS  NEGATIVE  COSTS 
C 

RMIN=0. 

DO  50  1=1, N 

IF(C(  D.GE.E)  GO  TO  50 

IF(C(  D.GE.RMIN)  GO  TO  50 

RMIN=C( I) 
50  CONTINUE 

IFIRMIN.EO.C. )  GO  TO  60 

DO  52  1=1, N 

52  C(I)=C(  D-RMIN 
N=N+1 

M  =  M+1 

H(M)=1000.*FF 
C(N)=-RMIN 
DO  53  J=1,N 

53  A(M,J)=FF 
C 

C      INITIALIZING  THE  TABLEAU 
C 

60  DO  3  J=1,M 
DO  2  1=1, M 

2  B(I,J)=0. 

3  B(J,J)=1. 
DO  4  1=1, M 
DB(I  )=-l. 
DUAL! I )=0. 

4  IB( I )=-I 
DO  5  J=1,N 
R(J)=C. 

D( J)=C. 

DO  5  1=1, M 

D( J)=D( J)-A( I, J) 

5  CONTINUE 

w=o» 

DO  6  1=1, M 

6  W=W+H(I) 
ITR=0 

C 

C      OBTAINING  A  COST  EQUAL  TO  ZERO  AND  A  PIVOT  COLUMN 

7  IF(DABS(W)-E)100,100,8 
3  DO  9  J=1,N 

IF(DU)  )10,9,9 
9  CONTINUE 
GO  TO  101 
10  JS=0 

DO  11  J=1,N 

IF(R(J).E0.1.)  GO  TO  11 
IF(DU).GE.-E)  GO  TO  11 
IF(JS.EO.O)  GO  TO  13 
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c 
c 
c 


12  IF(CX+C(J)/D( J) )11,11,13 

13  CX=-C(J)/D( J) 
JS  =  J 

11  CONTINUE 
47  DO  14  J=1,N 

14  C( J)=C( J)+CX*D( J) 
C( JS»=0. 

DO  301  J=1,N 
IF(CU).LT.E)  GO  TO  301 
GO  TO  302 

301  CONTINUE 
GO  TO  100 

302  Z=0. 

DO  303  I=lfM 

J=IB(I) 

IF(J.LE.C)  GO  TO  303 

Z=Z+CE( J)*H( I) 

303  CONTINUE 

DO  40  1=1, M 

40  DUAL( I)=DUAL( I )-CX*0B(  I) 

15  IF(W-E)100,100,41 

41  DO  16  J  =  1,N 

IF(C( J) -0.00001) 17,17,16 

17  C(J)=0. 

18  IF(0( J)-D< JS) )19,16,16 

19  JS=J 

16  CONTINUE 
WRITE  (6,213)  JS,  D(JS) 


213  FORMAT  (1H,  'JS  = 


20 
26 


100 
204 

101 
205 

425 
426 
526 


IF(DIJS) )2C,7,7 
CALL  ELEM 
IF(DI-1. )528,26,26 
CALL  PIVOT 
GO  TO  15 

PRINTING  OF  RESULTS 


14,  3X 


D= 


F12.4) 


429 

427 
428 
527 

528 


WRIT 
FORM 
GO  T 
WRIT 
FORM 
STOP 
DO  4 
WRIT 
FORM 
'VAL 
•DUA 
IFJR 
ZZ=0 
DO  4 
J=IB 
IF(J 
ZZ  =  Z 
CONT 
GO  T 
ZZ=Z 
WRIT 
FORM 
IF(I 
RETU 


E  (6,204) 

AT(1H  ,  'RESULT') 

0  425 

E  (6,205) 

AT( 1H  ,  • INFEASIBLE'  ) 

26  1=1, M 

E  (6,526)  IB(I),  H(I),  DUAL(I) 

AT  (1H  ,  'BASIS  VARIABLE  NUMBER=', 
UE  OF  BASIS  VARIABLES,  F12.4,  5X, 
L=',  F12.4) 
MIN.EO.O. )  GO  TO  427 

29  K=1,M 

(K) 

.LE.O)  GO  TO  429 

Z+CE< J)*H(K) 

INUE 

0  428 


14,  5X 


E  (6,527)  W,  ZZ 

AT  (1H  ,  'W=»,  F12.8, 

NDEX.EO.l)  GO  TO  1 

RN 


5X,  'Z=«,  F12.5) 
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c 

C      THE  NEW  TECHNIQUE 


C 


1  IT  =  0 

00  30  1=1, M 

IF(H(  D.GT.-E)  GO  TO  30 

IF(IT.EO.O)  GO  TO  31 

IF(H( I).GE.H( IT) )  GO  TO  30 
31  IT=I 
30  CONTINUE 

IF(IT.EO.O)  GO  TO  528 

IB(IT)=-IT 

H(IT)=-H( IT) 

DO  34  K  =  1,M 

DB(K)=B(  IT,K) 

34  B( IT,K)=-B(  IT,K) 
00  33  1=1, M 

IF(H(  D.GT.-E)  GO  TO  33 
H(I)=H( I)+H(  IT) 
DO  35  K=1,M 

35  B(I,K)  =  B(  I,K)+B(  IT,K) 
33  CONTINUE 

W=H( IT) 
DO  36  J=l,N 
AX(IT)=0. 
DO  21  K=1,M 
21  AX(IT)  =  AXUT)+B<IT,K)*A(Kf  J) 

36  D( J)=-AX( IT) 
INDEX=0 

GO  TO  8 
END 
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SUBROUTINE  ELEM 


IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON  4(70,150),B(70,70).C(15C) ,0(150) ,H(70) ,AX(70) 

COMMON  OB (70) ,DUAL(70),R( 1 50 ) ,CE ( 1 50 ) , E , W, Z , 0  I , FF , EP 

COMMON  CNAME( 150) 

COMMON  IB(7C),M,N,JS, I R , INDEX , I TR , Ml ,N1 

IR  =  0 

J  =  JS 
C 

C      CALCULATING  THE  REVISED  SIMPLEX  TABLEAU  PIVOT  COLUMM 
C      WITH  ROUND-OFF  ERROR  CONTROL 
C 

DO  22  1=1, M 

DM=0. 

DP=0. 

AX(I )  =  0. 

DO  21  K  =  1,M 

T=B( I,K)*A(K,J) 

IF(T.GE.O. )  GO  TO  40 

DM=DM+T 

GO  TO  21 
40  DP=DP*T 

21  CONTINUE 
AX(I )=DP+DM 
IF(DP.EQ.O. )  GO  TO  22 
IF(DABS(AX( I )/DP).GT.EP)  GO  TO  2  2 
AX(I  )=0. 

22  CONTINUE 
C 

C      FINDING  THE  PIVOT  ROW 


C 


JR=0 
EE  =  E 
31  DO  25  1=1, M 

GO  TO  5 


uu  ^5  l=i,M 

IF(INDEX.EO.O)  GO 

IF ( I B( I )  )  5,25,25 

IF(AX(I)-E)25,25,25 

IF(IR.EQ.O)  GC  TO  24 

IF(DABS(H( I )-AX( I  )  *X ) -EE ) 70, 70 , 2 

IF(H(  I)-AX(  I)*X)24,70,25 

X=H( I )/AX( I ) 


2 

24 
3  IR=I 

EE=E/AX( IR) 
GO  TO  25 
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c 
c 


SUBROUTINE  PIVOT 


IMPLICIT  REAL*8  (A-H,0-Z) 
COMMON  A (70, 15C  )  ,  B ( 70 , 70 ) ,C < 15C)tD( 150) ,H(70) , AX( 701 
COMMON  DB<70) .DUAL (70 >,R(150),CF(150),F,W,Z,DI,FF,EP 
COMMON  CNAME(150) 

COMMON  IB(70),M,N,JS,IR,INDFX,ITR,M1,N1 
OIMENSION  P(200) ,BX(200) 
WRITE  (6,214)  IR,  AX(IR) 
214  FORMAT  (1H,  25X,  MR=«,  13,  3X,  «AX=',  F12.M 

C      PERFORMANCE  OF  THE  PIVOT 
C 

43  IF(IBUR))  1,1,2 

2  K=IB(IR) 
R(K)=C. 
GO  TO  3 

1  P(JS)=1. 

3  IB(IR)=JS 
DI=DI*DABS( AX( IR) ) 
IFtDI.GE. 9999998.)  GO  TO  40 
DI  =  IDINT(DH-0.5) 

40  CONTINUE 
E=0.5/DI 
DO  27  1=1, M 

27  P(  I)=-AX( I ) /AX( IR) 
P( IR)=1./AX( IR)-1. 
P(M+1)=-D(JS)/AX( IR) 
DO  28  J=1,M 
BX(J)=B( IR,J) 

28  CONTINUE 
HX=H(  IR) 

00  29  J=1,M 

DO  29  1  =  1, M 

T=B( I, J) 

B(I, J)=B( I,J)+P(I)*BX( J) 

IF(T.EO.O. )  GO  TO  29 

IF(DABS(B(I,J)/T).GT.EP)  GO  TO  29 

B( I , J)=0. 

29  CONTINUE 

DO  30  1=1, M 

T=H( I) 

H(I)=H( I)+P< I)*HX 

IF(T.EO.O. )  GO  TO  3C 

IF(DABS(H(  I)/T).GT.EP)  GO  TO  30 

H ( I ) —0 • 

30  CONTINUE 
W=W-P(M4-1)*HX 
00  31  J=1,M 
T=DB( J) 

DB( J)=DB(J)+P(M+L)*BX(J) 

IF(T.EO.O. )  GO  TO  31 

IF(OABS(DB( J)/T).GT.EP)  GO  TO  31 

DB(J)=0. 

31  CONTINUE 
C 

C      UPDATING  THE  DB,  B,  AND  H  VALUES 

IF(DI. GE. 999998. )  GO  TO  81 
DO  92  J=1,M 

DF=DI*DABS(DB(J))*0.5 
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IF(DF.GT. 999998.)  GO  TO  61 
OF=IDINT(DF)/ni 
D8( J)=OSIGN(DF, D8( J) I 
61  00  92  1=1, M 

DF  =  DI*DABS(B(  I, J) )+0.5 

IF(DF. GT. 999998. )  GO  TO  92 

DF=IDINT(DF)/DI 

B(I, J)=DSIGN(OF,B( I, J) ) 

92  CONTINUE 

00  93  1=1, M 

DF=DI*DABS(H(  I  )  ) +0 . 5 

IF(DF.GT. 999998.)  GO  TO  93 

CF=IDINT(DF)/DI 

H( I)=DSIGN(DF,H(  I  )  ) 

93  CONTINUE 
81  CONTINUE 

DF=W*DI 

IF(DF. GE. 9999998. )  GO  TO  87 

W=IDINT(DI*W+0.5)/DI 
87  CONTINUE 
C 

C      CALCULATION  OF  0<J)»S  FOR  NEXT  ITERATION 
C 

DO  34  J=1,N 

DM=0. 

DP=0. 

D( J)=0. 

00  3  3  K  =  1,M 

T=DB(K)*A(K,J) 

IFU.GE.O.)  GO  TO  41 

OM=DM+T 

GO  TO  33 
41  DP=DP+T 

33  CONTINUE 
D( J)=DP+DM 

IF(DP.EO.O. )  GO  TO  34 
IF(DABS(D( J)/DP).LE.EP)  GO  TO  35 
IF(0A8S(D(J) )-E)35,35,34 

35  D(J)=0. 

34  CONTINUE 
D( JS)=0. 

IFCW.GT. 0.001)  GO  TO  44 
W=0 
44  ITR=ITR+1 

WRITE  (6,210)  ITR,  W,  Z,  01 
210  FORMAT  (1H  ,  5X,  MTR=»,  15,  5X,  • W= • ,  F18.8, 
15X,  «Z=»,  F12.4,  5X,  »0!*»,  E12.4) 

20  RETURN 

21  STOP 
FND 
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